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Abstract 

Modeling the flow of Non-Newtonian fluids in porous media is a challenging sub- 
ject. Several approaches have been proposed to tackle this problem. These include 
continuum models, numerical methods, and pore-scale network modeling. The lat- 
ter proved to be more successful and realistic than the rest. The reason is that it 
captures the essential features of the flow and porous media using modest computa- 
tional resources and viable modeling strategies. In this article we present pore-scale 
network modeling techniques for simulating non-Newtonian flow in porous media. 
These techniques are partially validated by theoretical analysis and comparison to 
experimental data. 
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1 Introduction 

The flow of non-Newtonian fluids in porous media is highly important subject and 
has many applications as these fluids are the rule rather than the exception. These 
applications include filtration of polymer solutions, enhanced recovery from oil 
reservoirs, medical and biological technologies, and soil remediation by removing 
toxic substances. Newtonian fluids are those fluids exhibiting a direct propor- 
tionality between stress and strain rate in laminar flow. All fluids for which this 
proportionality is violated, due to nonlinearity or initial yield-stress, are said to 
be non-Newtonian. These fluids are commonly divided into three broad groups: 
time-independent in which strain rate solely depends on the instantaneous stress, 
time-dependent in which strain rate is a function of both magnitude and duration 
of the applied stress, and viscoelastic which shows partial elastic recovery on re- 
moval of the deforming stress. A large number of models have been proposed in the 
literature to model the bulk rheology of non-Newtonian fluids under various flow 
conditions. A range of modeling strategies and computational techniques have also 
been developed to depict the flow of these fluids in porous media. In this article we 
highlight the computational techniques that have been developed and implemented 
in our non-Newtonian code using network modeling approach. 

In the context of fluid flow, ‘porous medium’ can be defined as a solid matrix 
through which small interconnected cavities occupying a measurable fraction of 
its volume are distributed. These cavities are of two types: large ones, called 
pores and throats, which contribute to the bulk flow of fluid; and small ones, 
comparable to the size of the molecules, which do not have an impact on the bulk 
flow though they may participate in other transportation phenomena like diffusion. 
The mathematical description of the flow in porous media is extremely complex 
and involves many approximations. In this regard, various methodologies have 
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been proposed and used. These include the continuum models, numerical methods 
and pore-scale network modeling. In this article we focus on network modeling as 
it is the methodology that we adopted in our flow simulation techniques. 

Pore-scale network modeling is a relatively novel method developed to deal with 
the flow through porous media and other related issues. It can be seen as a com- 
promise between the two extremes of continuum and numerical approaches as it 
partly accounts for the physics of flow and void space structure at pore level using 
affordable computational resources. Network modeling can be used to describe a 
wide range of properties from capillary pressure characteristics to interfacial area 
and mass transfer coefficients. The void space is described as a network of flow 
channels with idealized geometry. Rules that determine the transport properties 
in these channels are incorporated in the network to compute effective properties 
on a mesoscopic scale. The appropriate pore-scale physics combined with a repre- 
sentative description of the pore space gives models that can successfully predict 
average behavior [4, 5]. 

The general feature of network modeling is the representation of pore space 
by a network of interconnected ducts (bonds or throats) of regular shape and 
the use of a simplified form of the flow equations to describe the flow through the 
network. A numerical solver is normally employed to solve a system of simultaneous 
equations to determine the flow field. The network can be two-dimensional or 
three-dimensional with a random or regular lattice structure such as cubic. The 
shape of the cylindrical ducts include circular, square and triangular cross section 
and may include converging-diverging feature. The network elements may contain, 
beside the conducting ducts, nodes (pores) that can have zero or finite volume and 
may well serve a function in the flow phenomena or used as junctions to connect 
the bonds. The simulated flow can be Newtonian or non-Newtonian, single-phase, 
two-phase or even three-phase. The physical properties of the flow and porous 
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medium that can be obtained from flow simulation include absolute and relative 
permeability, formation factor, resistivity index, volumetric flow rate, apparent 
viscosity, threshold yield pressure and much more. Typical size of the network is a 
few millimeters. One reason for this minute size is to reduce the computational cost. 
Another reason is that this size is sufficient to represent a homogeneous medium 
having an average throat size of the most common porous materials. Up-scaling 
the size of a network is a trivial task if larger pore size is required. Moreover, 
extending the size of a network model by attaching identical copies of the same 
model in any direction or imposing repeated boundary conditions is another simple 
task. 

The general strategy in network modeling is to use the bulk rheology of the fluid 
and the void space description of the porous medium as an input to the model. The 
flow simulation in a network starts by modeling the flow in a single capillary. For a 
network of capillaries, a set of equations representing the capillaries and satisfying 
mass conservation have to be solved simultaneously to ford the pressure held and 
other physical quantities. For a network with n nodes there are n equations in n 
unknowns. These unknowns are the pressure values at the nodes. The essence of 
these equations is the continuity of how of incompressible fluid at each node in the 
absence of sources and sinks. To hnd the pressure held, this set of equations have 
to be solved subject to the boundary conditions which are the pressures at the inlet 
and outlet of the network. This unique solution is ‘consistent’ and ‘stable’ as it 
is the only mathematically acceptable solution to the problem, and, assuming the 
modeling process and the mathematical technicalities are reliable, should mimic the 
unique physical reality of the pressure held in the porous medium. For Newtonian 
fluid, a single iteration is needed to solve the pressure held as the how conductance 
is known in advance because the viscosity is constant. For purely viscous non- 
Newtonian fluid, the process starts with an initial guess for the viscosity, as it is 



1 INTRODUCTION 


unknown and pressure-dependent, followed by solving the pressure field iteratively 
and updating the viscosity after each iteration cycle until convergence is reached. 
For memory fluids, the dependence on time must be taken into account when 
solving the pressure field iteratively. 



2 MODELING THE FLOW IN PORO US MEDIA 


9 


2 Modeling the Flow in Porous Media 


In our model we use three-dimensional networks built from a topologically-equivalent 
three-dimensional voxel image of the pore space with the pore sizes, shapes and 
connectivity reflecting the real medium. Pores and throats are modeled as having 
triangular, square or circular cross-section by assigning a shape factor which is 
the ratio of the area to the perimeter squared and obtained from the pore space 
image. Most of the network elements are not circular. To account for the non- 
circularity when calculating the volumetric flow rate analytically or numerically 
for a cylindrical capillary, an equivalent radius R eq is defined 


where G is the geometric conductance which may be obtained empirically from nu- 
merical simulation. The network can be extracted from voxel images of real porous 
materials or from voxel images generated by simulating the geological processes 
by which the porous medium was formed. Examples for the latter are the two 
networks of Statoil which represent two different porous media: a sand pack and 
a Berea sandstone. These networks are constructed by 0ren and coworkers [6, 7] 
and have been used by several researchers in flow simulation studies. Another pos- 
sibility for generating a network is by employing computational algorithms based 
on numeric statistical data extracted from the porous medium of interest. Other 
possibilities can also be found in the literature. An important aspect that charac- 
terizes the flow in porous media and makes it distinct from bulk is the presence 
of converging-diverging flow paths. This geometric factor significantly affects the 
flow and accentuates elastic responses. Therefore, a converging-diverging feature 
is introduced to the network capillaries when modeling viscoelastic flow. 



( 1 ) 


Assuming a laminar, isothermal and incompressible flow at low Reynolds num- 
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ber, the only equations that require attention are the constitutive equation for the 
particular fluid and the conservation of volume as an expression for the conserva- 
tion of mass. For Newtonian flow, the pressure field can be solved once and for all. 
For non-Newtonian flow, the situation is more complex as it involves non-linearities 
and requires iterative techniques. For the simplest case of time-independent fluids, 
the strategy is to start with an arbitrary guess. Because initially the pressure drop 
across each network element is not known, an iterative method is used. This starts 
by assigning an effective viscosity to the fluid in each element. The effective vis- 
cosity is defined as that viscosity which makes Poiseuillc’s equation fits any set of 
laminar flow conditions for time-independent fluids [8]. By invoking the conserva- 
tion of volume for incompressible fluid, the pressure field across the entire network 
is solved using a numerical solver [9]. Knowing the pressure drops in the network, 
the effective viscosity of the fluid in each element is updated using the expression 
for the flow rate in a capillary with the Poiseuillc’s law as a definition. The pressure 
field is then recomputed using the updated viscosities and the iteration continues 
until convergence is achieved when a specified error tolerance in the total flow rate 
between two consecutive iteration cycles is reached. Finally, the total volumetric 
flow rate and the apparent viscosity, defined as the viscosity calculated from the 
Darcy’s law, are obtained. Other physical parameters of interest that characterize 
the fluid and the porous medium may also be computed in this process. 
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3 Modeling Time-Independent Flow 

Three time-independent fluids have been incorporated in the non-Newtonian model. 
These are Carreau, Ellis and Herschel-Bulkley. Carreau is a four-parameter model 
given by 


where /i is the fluid viscosity, ^ is the viscosity at infinite shear, /i 0 is the viscosity 
at zero shear, 7 is the shear rate, n is the flow behavior index, and 7 cr is a critical 
shear rate given by 


where C is the consistency factor of the equivalent shear-thinning fluid in the power- 
law formulation. This model was previously implemented and fully described by 
Lopez [10] and Lopez et al [11]. In summary, the implementation of Carreau model 
relies on the use of an empirical expression for the volumetric flow rate in a single 
tube. The reader should refer to those references for details. 

Ellis is a three-parameter model which describes time- independent shear-thinning 
yield-free non-Newtonian fluids. It is used as a substitute for the power-law and 
is appreciably better than the power-law model in matching experimental mea- 
surements. Its distinctive feature is the low-shear Newtonian plateau without a 
high-shear plateau. According to this model, the fluid viscosity /i is given by [12- 
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where /i 0 is the low-shear viscosity, r is the shear stress, r 1/2 is the shear stress at 
which fi = Ho / 2 and a is a dimensionless indicial parameter related to the slope 
in the power-law region. For Ellis fluids, the volumetric flow rate in a circular 
cylindrical tube is given by [12-15]: 


nR 4 AP 4 / RAPY 1 

8 Lfi 0 a + 3 Lt 1/2 J 


( 5 ) 


where R is the tube radius, A P is the pressure drop across the tube and L is the 
tube length. 

Herschel-Bulkley is a three-parameter model that can describe Newtonian and 
a large group of time-independent non-Newtonian fluids. It is given by [8] 


t = r 0 + Cy n (r > T 0 ) 


( 6 ) 


where r is the shear stress, r D is the yield-stress above which the substance starts 
flowing, C is the consistency factor, 7 is the shear rate and n is the flow behavior 
index. Herschel-Bulkley reduces to the power-law, or Ostwald-de Waele model, 
when the yield-stress is zero, to the Bingham plastic model when the flow behavior 
index is unity, and to the Newton’s law for viscous fluids when both these conditions 
are met. For Herschel-Bulkley fluids, the volumetric flow rate in a cylindrical 
capillary at yield is given by [8]: 

n = — (—X(t — T 1 1+ " V Tw ~ To "> 2 I 27 °(^-7q) T q 2 
J c s \&P ) W 0 [ 3 + l/n + 2 + 1/n + 1 + l/n 

(' T w > T 0 ) 


( 7 ) 
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where t w (= is the shear stress at the tube wall. 

The implementation of Ellis and Herschel-Bulkley is based on the use of the 
analytical expressions for the volumetric flow rate in a circular cylindrical duct, as 
given above (i.e. Equation 5 for Ellis and Equation 7 for Herschel-Bulkley). These 
expressions are used in conjunction with an iterative technique to find the total 
flow across the network, as outlined in § 2. 
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4 Modeling Yield-Stress 


Yield-stress or viscoplastic fluids can sustain shear stresses, that is a certain amount 
of stress must be exceeded before the flow starts. So an ideal yield-stress fluid is 
a solid before yield and a fluid after. Accordingly, the viscosity of the substance 
changes from an infinite to a finite value. However, the physical situation suggests 
that it is more realistic to regard a yield-stress substance as a fluid whose viscos- 
ity as a function of applied stress has a discontinuity as it drops sharply from a 
very high value on exceeding a critical stress. Several constitutive equations to 
describe yield-stress substances are in use; the most popular ones are Bingham, 
Casson and Herschel-Bulkley. In our network model, yield-stress was implemented 
within the Herschel-Bulkley fluid. A number of numerical algorithms, related to or 
independent of Herschel-Bulkley, were also implemented in the model. 

The implementation of the yield-stress in a network is based on the yield con- 
dition for its conducting ducts which are assumed to be circular cylinders. The 
verification of the yield condition in the individual ducts associates the process 
of solving the pressure field in the network. For yield-stress fluids, the threshold 
pressure drop above which the flow in a single tube starts is given by 


A P th = 


2 Lt q 
R 


( 8 ) 


where A P th is the threshold pressure drop, r G is the yield-stress and R and L are 
the tube radius and length respectively. 

In our model, the substance before yield is considered to be fluid with very high 
but finite viscosity so the flow virtually vanishes. The reason is that the pressure 
across the network have to communicate. Accordingly, the pressure field in the case 
of yield-stress fluids is solved as in the case of non-yield-stress fluids. A further 
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condition is also imposed before any duct is allowed to yield, that is the duct must 
be part of a non-blocked path spanning the network from the inlet to the outlet. 
The logic is that any conducting duct should have a source on one side and a sink 
on the other. 

4.1 Predicting Threshold Yield Pressure of a Network 

Non-Newtonian literature contains several attempts to predict the yield point of 
a complex porous medium from the void space description and yield-stress value 
of an ideal yield-stress fluid without modeling the flow process. In this regard, 
there is an implicit assumption that the network is an exact replica of the medium 
and the yield-stress value reflects the real yield-stress of fluid so that any failure 
of these proposals can not be attributed to mismatch or any factor other than 
flaws in these proposals. Our discussion in this section will focus on the Invasion 
Percolation with Memory (IPM) and Path of Minimum Pressure (PMP) algorithms 
which are implemented in the network model to make such predictions. 

The IPM is an algorithm for finding the inlet-to-outlet path that minimizes the 
sum of values of a property assigned to the individual elements of the network, 
and hence finding this minimum. For a yield-stress fluid, this reduces to finding 
the inlet-to-outlet path that minimizes the yield pressure. The yield pressure of 
this path is then taken as the network threshold yield pressure. A flow chart of 
the IPM algorithm is presented in Figure (1). The PMP algorithm is based on 
a similar assumption to that upon which the IPM is based, that is the network 
threshold yield pressure is the minimum sum of the threshold yield pressures of 
the individual elements of all possible paths from the inlet to the outlet. However, 
PMP is computationally different and is more efficient than the IPM in terms 
of the required computational resources (memory and CPU time). A flow chart 
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depicting the PMP algorithm is given in Figure (2). In most cases that have been 
investigated the IPM and PMP produce identical results. 

Another algorithm, called Actual Threshold Pressure (ATP), for finding the 
network yield point by solving the pressure field was also developed and imple- 
mented. The ATP is an iterative simulation algorithm which uses the solver to 
find the network yield pressure to the required accuracy. However, this algorithm 
is not independent of Herschel-Bulkley model as it relies on multiple applications 
of flow simulation of this model. 
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Figure 1: Flowchart of the Invasion Percolation with Memory (IPM) algorithm. 
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Figure 2: Flowchart of the Path of Minimum Pressure (PMP) algorithm. 
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5 Modeling Time-Dependent Flow 

In theory, time-dependence effects can arise from thixotropic (and rheopectic) 
structural change or from time-dependent viscoelasticity or from both effects simul- 
taneously. The existence of these two different types of time-dependent rheological 
behavior is generally recognized. Although it is convenient to distinguish between 
these as two separate phenomena, real fluids can exhibit both types simultaneously. 
Several physical distinctions between viscoelastic and thixotropic time-dependence 
have been made. The important one is that while time-dependence of viscoelastic 
fluids arises from relaxation and delayed response, time-dependence of thixotropic 
fluids arises from structural change. With regards to modeling the flow in porous 
media of complex fluids that have time dependency in a dynamic sense due to 
thixotropic or clastic nature, there are three major difficulties 

• The difficulty of tracking the fluid elements in the pores and throats and 
identifying their deformation history, as the path followed by these elements 
is random and can have several unpredictable outcomes. 

• The mixing of fluid elements with various deformation history in the individ- 
ual pores and throats. As a result, the viscosity is not a well-defined property 
of the fluid in the pores and throats. 

• The change of viscosity along the streamline since the deformation history is 
continually changing over the path of each fluid element. 

These complications have not been considered in the current model, and hence 
no dynamic time-dependence has been included in the code. However, general 
strategies for simulating time-dependent thixotropic behavior have been considered. 
These strategies can provide a framework for future development. There are three 
major cases of flow simulation of thixotropic fluids in porous media: 
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• The flow of strongly strain-dependent fluid in a porous medium which is not 
highly homogeneous. This case is very difficult to model clue to the difficulty 
of tracking the fluid elements in the pores and throats and determining their 
deformation history. Moreover, the viscosity function may not be well defined 
due to the mixing of fluid elements with various deformation history in the 
individual pores and throats. 

• The flow of strain-independent or weakly strain-dependent fluid through 
porous media in general. A possible strategy is to apply single time-dependent 
viscosity function to all pores and throats at each instant of time and hence 
simulating time development as a sequence of Newtonian states. 

• The flow of strongly strain-dependent fluid in a highly homogeneous porous 
medium such that the fluid is subject to the same deformation in all ducts. 
The strategy for modeling this flow is to define an effective pore strain rate. 
Then using a very small time step the strain rate in the next instant of time 
can be found assuming constant strain rate. As the change in the strain rate 
is then known, a correction to the viscosity due to strain-dependency can be 
introduced. 

It should be remarked that the Bautista-Manero fluid, which is used to model 
steady-state viscoelastic flow, incorporates thixotropic as well as viscoelastic at- 


tributes. 
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6 Modeling Viscoelastic Flow 

As indicated already, no dynamic time-dependence has been included in the non- 
Newtonian flow model. For the steady-state flow of viscoelastic fluids, the approach 
of Tardy [16, 17] was used with some adaptation. In this approach, the capillaries 
of the network are modeled with contraction to account for the effect of converging- 
diverging geometry on the flow. The reason is that the effects of fluid memory take 
place on going through a radius change, as this change induces a change in strain 
rate with viscosity changing consequences. The capillaries are also discretized in 
the flow direction and a discretized form of the flow equations is used with assumed 
prior knowledge of the stress and viscosity at the inlet of the network. Starting with 
an initial guess for the flow rate and using iterative technique, the pressure drop as 
a function of the flow rate is then found in each capillary. Finally, the pressure field 
for the whole network is found iteratively until convergence is achieved. Once this 
happens, the flow rate through each capillary in the network can be computed and 
the total flow rate through the network is determined by summing and averaging 
the flow through the inlet and outlet capillaries. 

This algorithm employs a one- dimensional version of the Bautista-Manero model 
which combines the Fredrickson kinetic equation for flow-induced structural changes 
with the Oldroyd-B constitutive equation for viscoelasticity. The model requires six 
parameters that have physical significance and can be estimated from rheological 
measurements. Bautista-Manero model was originally proposed for the rheology of 
worm-like micellar solutions which usually have an upper Newtonian plateau, and 
show strong signs of shear-thinning. The model, which incorporates shear-thinning, 
elasticity and thixotropy, can be used to describe the complex rheological behavior 
of viscoelastic systems that also exhibit thixotropy and rheopexy under shear flow 
[16, 18-20]. 
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The kinetic equation of Fredrickson that accounts for the destruction and con- 
struction of structure is given by 


where /i is the viscosity, t is the time of deformation, A is the relaxation time, fi 0 and 
/i oo are the viscosities at zero and infinite shear rates respectively, k is a parameter 
that is related to a critical stress value below which the material exhibits primary 
creep, r is the stress tensor and 7 is the rate of strain tensor. In this model, A 
is a structural relaxation time, whereas k is a kinetic constant for structure break 
down [16, 18-20], 

The Oldroyd-B constitutive equation is given by [14] 


where Ai is the relaxation time, A 2 is the retardation time, and 7 is the upper 
convected time derivative of the rate-of-strain tensor given by 


Similar expression applies to the upper convected time derivative of the stress tensor 





( 10 ) 


V 



(ii) 


t. A flow chart outlining the steady-state viscoelastic flow simulation algorithm is 


presented in Figure (3). 
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Figure 3: Flow chart of the steady-state viscoelastic flow simulation algorithm. 
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7 Validation and Results 

These computational techniques as implemented in the non-Newtonian code were 
quantitatively validated in a number of cases. These include low flow rate regimes, 
Bingham fluids at high flow rates, Boger fluids, and Newtonian as a special case 
of non-Newtonian model [17, 21, 22], Two randomly-distributed networks repre- 
senting two different porous media, a sand pack and a Berea sandstone, were used 
in these validations. These networks are constructed by 0ren and coworkers [6, 7] 
from voxel images generated by simulating the geological processes by which the 
porous medium was formed. The sand pack comprises 13490 elements (pores and 
throats) with a cube side length of 2.5mm, while the Berea consists of 38495 ele- 
ments with a cube side of 3mm. The physical and statistical properties of these 
networks with detailed comparison between them can be found in [17, 21]. 

These computational techniques were also partially validated by a number of 
experimental data sets found in the literature for time-independent fluids [17, 21]. 
A sample of these data with their simulation counterparts are given in Figures (4), 
(5) and (6) for Ellis, Herschel-Bulkley and Bingham fluids respectively. The bulk 
rheologies of these data sets are presented in Tables (1), (2) and (3) respectively. 
In these simulations, scaled versions of the sand pack network were used. The 
purpose of scaling is to match the network characteristics to the characteristics of 
the corresponding porous media. The sand pack was used instead of Berea because 
it is a better match to the experimental packed beds in terms of homogeneity and 
tortuosity. The experimental validation is based on using the experimental bulk 
rheology and bed properties as an input to the non-Newtonian code. Qualitatively, 
all trends of behavior that have been observed are sensible. The major failure of 
the non-Newtonian model occurs in the case of fluids with yield-stress. Although 
the quality of some experimental data sets may be questionable, it seems that the 
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yield-stress model, which is based on the concept of equivalent radius of cylindrical 
capillaries, is too simplistic and unrealistic and hence is very unlikely to produce 
reliable predictions. 



Figure 4: Park’s Ellis experimental data sets [1] for polyacrylamide solutions with 
0.50%, 0.25%, 0.10% and 0.05% weight concentration flowing through a coarse 
packed bed of glass beads having K = 3413 Darcy and </> = 0.42 alongside the 
simulation results obtained with a scaled sand pack network having the same K 
presented on a log- log scale. The bulk rheology of these data is given in Table 1. 
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Figure 5: Al-Fariss and Pinder’s Herschel-Bulkley experimental data sets [2] for 
5.0% wax in Clams B oil flowing through a column of sand having K = 315 Darcy 
and (j) = 0.36 alongside the simulation results obtained with a scaled sand pack 
network having the same K and </>. The temperatures, T, are in °C. The bulk 
rheology of these data is given in Table 2. 



Figure 6: Network simulation results with the corresponding experimental data 
points of Chase and Dachavijit [3] for Bingham aqueous solutions of Carbopol 
941 with various concentrations (0.37%, 0.45%, 0.60%, 1.00% and 1.30%) flowing 
through a packed column of glass beads. The bulk rheology of these data is given 
in Table 3. 
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Table 1: The bulk rheology of Park’s Ellis experimental data sets [1]. 


Dataset 

ho (Pa.s) 

a 

D/2 ( Pa ) 

0.50% 

4.35213 

2.4712 

0.7185 

0.25% 

1.87862 

2.4367 

0.5310 

0.10% 

0.60870 

2.3481 

0.3920 

0.05% 

0.26026 

2.1902 

0.3390 


Table 2: The bulk rheology of Al-Fariss and Pinder’s Herschel-Bulkley experimental 
data sets [2], 


Dataset 

C (Pa.s n ) 

n 

t 0 (Pa) 

T=16 

0.463 

0.87 

3.575 

T=18 

0.568 

0.80 

2.650 

T=20 

0.302 

0.90 

1.921 


Table 3: The bulk rheology of Chase and Dachavijit’s Bingham experimental data 
sets [3]. 


Dataset 

C (Pa.s 11 ) 

n 

T 0 (Pa) 

0.37% 

0.017 

1.0 

2.06 

0.45% 

0.038 

1.0 

4.41 

0.60% 

0.057 

1.0 

7.09 

1.00% 

0.128 

1.0 

17.33 

1.30% 

0.215 

1.0 

28.46 
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The computer code in which these computational techniques have been imple- 
mented is derived from the non-Newtonian code of Lopez [10] which was con- 
structed from an early version of the Newtonian code by Valvatne [23]. Beside 
the Carreau model, which is inherited from the original code of Lopez, the current 
code can simulate the flow of Ellis and Herschel-Bulkley fluids. Several algorithms 
related to yield-stress and a modified version of the Tardy algorithm to simu- 
late steady-state viscoelastic flow using the Bautista-Manero model are also imple- 
mented. The code can be downloaded from this URL: http : //www3 . imperial . ac . 
uk/earthscienceandengineering/research/perm/porescalemodelling/ software/ 
non-newtonian°/ 0 20code or this URL: www.scienceware.net/id2.html. 

The code has a command line interface that uses a keyword-based input hie. 
Convergence time generally depends on the fluid rheology, the size of network and 
the type of algorithm. A typical convergence time for the sand pack and Berea 
sandstone networks used in this study is a second for the time- independent models 
and a few seconds for the viscoelastic model. The time requirement for the yield- 
stress algorithms is highly dependent on the last two factors. However, there is 
significant difference between the IPM and PMP convergence time. As these two 
algorithms produce very similar results, it is recommended to use the PMP for 
large networks. In all cases, the memory requirement does not exceed a few tens 
of megabytes for a network with up to 12000 pores. In general, the memory cost is 
affordable on a typical modern workstation for all available networks. The general 
how sequence of the program is as follows: 

• The program starts by reading the input and network data hies followed by 
creating the network. 

• For fluids with yield-stress, the program executes IPM and PMP algorithms 
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to predict the threshold yield pressure of the network. This is followed by 
an iterative simulation algorithm to find the network actual threshold yield 
pressure to the required accuracy. 

• Single-phase flow of a Newtonian fluid is simulated to find the fluid-related 
network properties such as absolute permeability. 

• Single-phase flow of Newtonian and non-Newtonian fluids is simulated for a 
range of pressure points as defined by the user. 

In all stages, informative messages are issued about the program progress and 
the data obtained. The program also creates several output data files. These 
include script files to visualize the entire network or slice of it and the flow path of 
yield-stress fluid using Rhino 3D program. 
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9 Conclusions and Discussions 

In this study, we outlined a set of computational techniques based on pore-scale 
network modeling to simulate single-phase flow of non-Newtonian fluids in porous 
media. These techniques are implemented in a computer code and have been 
partially validated analytically and experimentally. 

1. The success was evident in the case of time-independent fluids. This includes 
comparison with a number of experimental data sets and correct predictions 
in special and limiting cases such as Newtonian fluids and the asymptotic 
behavior of Bingham at high flow rates. 

2. Steady-state viscoelastic flow simulation has also produced reliable results in 
the case of low flow rate regimes and Boger fluids. Moreover, qualitatively 
sensible trends of behavior were observed in the other cases using several 
parameters related to the fluid, porous media and numerical indicators. 

3. Thixotropic flow as such has not been modeled although thixotropic aspects 
are included within the Bautista-Manero model which is the basis of the 
steady-state viscoelastic flow algorithm. However, thixotropic computational 
strategies have been developed and assessed. 

4. The predictions were less satisfactory in the case of yield-stress fluids. This 
may be explained by inadequate representation of the pore space structure, 
experimental errors and involvement of other physical phenomena. Minimum 
threshold path algorithms (i.e. IPM and PMP) have also been developed and 
implemented. The analysis revealed that these algorithms are too simplistic 
and hence cannot produce reliable predictions for the pressure yield point of 


a network. 
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Nomenclature 


a 

parameter in Ellis model 


7 

strain rate (s _1 ) 


7er 

critical shear rate (s _1 ) 


7 

rate-of-strain tensor 


A 

structural relaxation time in Fredrickson model (s) 


Ai 

relaxation time (s) 


A2 

retardation time (s) 


9 

viscosity (Pa.s) 


9 0 

zero-shear viscosity (Pa.s) 


Loo 

infinite-shear viscosity (Pa.s) 


T 

stress (Pa) 


T 

stress tensor 


T i/ 2 

stress when /j = /j 0 /2 in Ellis model (Pa) 


To 

yield-stress (Pa) 


Tw 

stress at tube wall (Pa) 



porosity 


c 

consistency factor (Pa.s n ) 


G 

geometric conductance (m 4 ) 


G' 

flow conductance (m 3 .Pa _1 .s _1 ) 


k 

parameter in Fredrickson model (Pa^ 1 ) 


K 

absolute permeability 


L 

tube length (m) 
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n flow behavior index 

P pressure (Pa) 

A P pressure drop (Pa) 

A P t h threshold pressure drop (Pa) 

Q volumetric flow rate (m 3 .s -1 ) 

R tube radius (m) 

R eq equivalent radius (m) 

t time (s) 

T temperature (K, °C) 

v fluid velocity vector 

ATP Actual Threshold Pressure algorithm 

IPM Invasion Percolation with Memory algorithm 

PMP Path of Minimum Pressure algorithm 

upper convected time derivative 
(•) T matrix transpose 


Note: units, when relevant, are given in the SI system. Vectors and tensors are 
marked with boldface. Some symbols may rely on the context for unambiguous 


identification. 
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